碎碎念一下

是說，由於第一天想起 #鐵人賽 開賽這件事已經是晚上11點40幾分的事情了，待我打完第一天的文章時剛好已過12點，結果就變成10/2的發文了 QQ

安裝好開發Julia code的環境之後

那我們就趕緊先進入Julia，順便安裝一下幾個我馬上就會用到的packages： DataFrames、CSV跟BioSequences。這邊我們就先不管基礎語法，先管安裝我們需要用到的東西即可。

在Julia裡面，當我們需要安裝、更新移除packages時，我們絕對會使用到的一個Package叫做Pkg，在引入了這個package之後，我們就可以它提供的各種函數來操作各個Package。

using Pkg
pkgs = ["DataFrames", "CSV", "BioSequences"]
Pkg.add(pkgs)

沒什麼意外的話，藉由上面的操作，稍微等待一些時間，就可以順利完成這三個Packages的安裝。在今天的牛刀小試中，我打算拿Julia來處理這兩天我正在處理的SNP array [1] data。雖然這幾年NGS技術進展地很快，定序的價格也是直直落（見下圖），但在臨床端或是一些比較大規模的基因體變異研究當中，為了節省成本，大部分還是會選擇使用晶片來作為偵測人體基因體變異位點的主要工具。

這週開始正好得開始處理公司所接到的一個合作案，目標也是想要解決一些population genetics上的問題。由於先前完全沒有相關的經驗，在爬了幾篇文獻跟網路上的介紹文章之後，注意到了有不少研究都會提到，由於那些SNP array平台商所提供的自動化軟體裡面用來做genotyping (基因分型)的演算法不完全正確，且這些晶片在掃描螢光訊號的時候，疑似(?)也很容易發生一些把訊號強度很接近的probes分到同一群裡面，造成最後的genotyping的結果會有錯誤，因此這些文獻都不約而同地提到了要對raw intensity data進行re-calling一次，然後去掉那些不可信的probes。要進行re-calling，首先就得拿到所使用的那個晶片平台所有probes的序列，然後再將這些序列一一比對回human genome sequence。這次手上拿到的data是使用Illumina這間公司所提供的Human CoreExome-24 v1晶片平台做出的，為了提取出所有的探針序列，我們必須先從官網下載它們所提供的annotation資訊，我下載的是這個---HumanCoreExome-24v1-0_A.csv。

拿到這個檔案後，我們得先對這個檔案做一些處理，並將結果存在annotation.csv這個檔案中。這邊我習慣使用awk來對biological data做些前處理，相信使用linux的人也都很熟悉：

awk -F "," 'NR > 7 {print $0}' HumanCoreExome-24v1-0_A.csv | grep -v ^00 | grep -v Controls > annotation.csv

接下來，我得parse一下這個新產生的annotation.csv檔，一是要產生給各probe使用的ID，二是要提取出裡面的probe序列。另一方面，由於不是所有的probes都提供了兩種序列，所以必須分開提取。這邊，第一步我使用了CSV.read來讀取剛剛處理產生的annotation.csv，接下來根據AlleleB_ProbeSeq這個欄位有沒有值來判斷這個probe究竟是只有一種序列還是有兩種序列，而這個判斷會需要使用到ismissing這個function。最終，我把提取出來的序列使用BioJulia計畫中的BioSequences.jl這個package來存成FASTA格式的文件輸出。

#!/usr/bin/env julia
## Importing pkgs
using DataFrames
using CSV
using BioSequences

##  Inputs
annotFile = "annotation.csv"
annotation = CSV.read(annotFile, header=1, delim=',')

## Generating Fasta identifier
fastaID = [join(Matrix(annotation[i,1:2]), "_SEQID_") for i in 1:nrow(annotation)]

## Check probe status
probeA_Seq = annotation[6]
probeB_Seq = annotation[8]
nothasProbeB = ismissing.(probeB_Seq)
hasBothProbes = [!b for b in nothasProbeB]

# Select data without probeB info 
onlyProbeA_FastaID = fastaID[nothasProbeB]
onlyProbeA_Seq = probeA_Seq[nothasProbeB]
onlyProbeA_output = "only_probeA.fasta"

# Select data with Probe A and Probe B both
withBothProbes_FastaID = fastaID[hasBothProbes]
withBothProbes_ProbeA_Seq = probeA_Seq[hasBothProbes]
withBothProbes_ProbeB_Seq = probeB_Seq[hasBothProbes]
withBothProbes_output = "with_both_probes.fasta"

# Output fasta
fh1 = FASTA.Writer(open(onlyProbeA_output, "w"))
for i in 1:length(onlyProbeA_FastaID)
    rec = FASTA.Record("ONLY_PROBE_A", onlyProbeA_FastaID[i], DNASequence(onlyProbeA_Seq[i]))
    write(fh1, rec)
end
close(fh1)

fh2 = FASTA.Writer(open(withBothProbes_output, "w"))
for i in 1:length(withBothProbes_FastaID)
    probeA_rec = FASTA.Record("PROBE_A", withBothProbes_FastaID[i], DNASequence(withBothProbes_ProbeA_Seq[i]))
    probeB_rec = FASTA.Record("PROBE_B", withBothProbes_FastaID[i], DNASequence(withBothProbes_ProbeB_Seq[i]))
    write(fh2, probeA_rec)
    write(fh2, probeB_rec)
end
close(fh2)

執行了這支程式碼，就可以看到有兩個FASTA檔產生了。今天先到這，明天再來細講DataFrame跟CSV這兩個Packages的使用。

Reference:

1. SNP簡介:http://bioinfo.cs.ccu.edu.tw/wiki/doku.php?id=single_nucleotide_polymorphism_snp